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H 1 ABSTRACT 

§: 

We examine the non-linear stability of the Wisdom-Holman (WH) symplec- 
^ ■ tic mapping applied to the integration of perturbed, highly eccentric (e ^ 0.9) 

two-body orbits. We find that the method is unstable and introduces artificial 
chaos into the computed trajectories for this class of problems, unless the step 
size is chosen small enough to always resolve periapse, in which case the method 
is generically stable. This 'radial orbit instability' persists even for weakly per- 
turbed systems. Using the Stark problem as a fiducial test case, we investigate 
the dynamical origin of this instability and show that the numerical chaos results 
from the overlap of step size resonances (cf. Wisdom & Holman 1992); interest- 
ingly, for the Stark problem many of these resonances appear to be absolutely 
q . stable. 

£j , We similarly examine the robustness of several alternative integration meth- 

ods: a regularized version of the WH mapping suggested by Mikkola (1997); the 

£> ! potential-splitting (PS) method of Lee et al. (1997); and two methods incor- 

porating approximations based on Stark motion instead of Kepler motion (cf. 
Newman et al. 1997| ). The two fixed point problem and a related, more general 



problem are used to comparatively test the various methods for several types of 
motion. Among the tested algorithms, the regularized WH mapping is clearly 
the most efficient and stable method of integrating eccentric, nearly-Keplerian 
orbits in the absence of close encounters. For test particles subject to both high 
eccentricities and very close encounters, we find an enhanced version of the PS 
method — incorporating time regularization, force-center switching, and an im- 
proved kernel function — to be both economical and highly versatile. We conclude 
that Stark-based methods are of marginal utility in A-body type integrations. 
Additional implications for the symplectic integration of A-body systems are 
discussed. 



Subject headings: chaos — methods: numerical — celestial mechanics 
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1. Introduction 

Symplectic integration schemes have become increasingly popular tools for the numerical 
study of dynamical systems, a result of their often high efficiency as well as their typical long- 
term stability (see, e.g., [Marsdcu, Patrick, fc Shadwick 1996 and the many references within). 



The Wisdom-Holman (WH) symplectic mapping in particular ( |Wisdom fc Holman 19911 ; cf. 
|Kinoshita, Yoshida, fc Nakai 199~1~|) has been widely used in the context of Solar System 
dynamics. However, the fact that this and other symplectic methods are, by construction, 
"finely tuned" can make them susceptible to performance-degrading ailments (much as high- 
order methods offer little benefit if the motion is not sufficiently smooth), and the stability 
of these methods for arbitrary systems and initial conditions is not completely understood. 
It would be prudent, therefore, to exercise caution when applying such schemes to systems 
entering previously unexplored dynamical states, and to ensure that adequate preliminary 
testing is undertaken regardless of the method's stability in previously considered problems. 
Recent galactic dynamics simulations by Rauch & Ingalls (1997), for example — which used 
the WH mapping — uncovered evidence of an instability in the method when applied to a 
particular class of problems: the integration of highly elliptical, nearly-Keplerian orbits in 
which the timestep is taken small enough to smoothly resolve the perturbation forces, but 
not so tiny as to explicitly resolve pericenter. Since particle motion in these simulations 
was extremely close to Keplerian near pericenter, and since the mapping itself is exact for 
Keplerian motion, there is no a priori reason why the method should have performed as 
poorly and unstably as was found. 

Recently several variations of the WH mapping have been proposed which aim to extend 
the range of applicability of the original method. The regularized WH mappings investigated 
by Mikkola (1997), for instance, appear promising in the context of elliptical motion. The 
potential-splitting (PS) method of Lee, Duncan, & Levison (1997) allows symplectic inte- 
gration of close encounters between massive bodies by adding a multiple-timestep algorithm 
similar to that of Skeel & Biesiadecki (1994). Unfortunately both of these schemes have 
limitations of their own; the former is unable to resolve close encounters, while the latter 
approach (like the original mapping) appears to be unstable when orbits are eccentric (cf. 
Puncan, Levison, fc Lee 1997|) . 



In this paper, we use a series of test problems based on perturbed two-body motion 
to analyze the stability of the WH mapping and several of its variants. In particular, we 
examine the reliability of the methods for test particles whose motion is either highly eccentric 
or subject to close encounters with the perturbers (or both). The plan of the paper is as 
follows. In the following section, the performance of the WH mapping at high eccentricities 



is investigated using the Stark problem (e.g., Pankowicz 1994| [Kirchgraber 197i|) — for which 
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the range of orbital eccentricities is easily controlled, and no close encounters occur — as the 
fiducial test case. The instability found in the integrated motion is then explained using 
complementary analytic and geometric arguments. In § |^, modified forms of the original 
mapping are described; similarly, in § ^ integrators based on Stark motion instead of Kepler 
motion are considered. Section |] uses the two fixed point problem (e.g., [Pars 19651) as well 
as a more general test problem (drawn from the area of galactic dynamics) to conduct a 
comparative performance analysis of the various algorithms. Both the Stark and two fixed 
point problems are fully integrable and analytically soluble in terms of elliptic functions and 
integrals, allowing a detailed assessment of the accuracy of the numerical results to be made. 
Concluding discussion is given in § ||. 



2. The Stark Problem 

2.1. Orbital Motion 

The Stark problem represents the motion of a test particle about a fixed Newtonian 
force center (i.e., the Kepler problem) subject in addition to a uniform force of constant 
magnitude and direction. The Hamiltonian (per unit mass) for the problem is given by 

2 r 

where x is the Cartesian position, p is its conjugate momentum (in this case, the physical 
velocity), r = |je|, M is the mass of the central object, and the constant vector S (the 'Stark 
vector') embodies the uniform field. Physical analogies include: the (classical) orbit of an 
electron about a fixed nucleus immersed in a uniform electric field; the trajectory about its 
parent of an artificial satellite with a continuously firing thruster; and the motion of a dust 
particle around a comet nucleus taking into account the radiation pressure of sunlight (cf. 
[Hamilton 1992| ; [Mignard 1982] ). There are three conserved quantities (and hence the motion 



is regular): the energy, E, the angular momentum component along the Stark vector, L ■ S, 
and a third integral, a (say), which arises as a 'separation constant' in the analytic solution 
of the problem. 

Qualitatively speaking, motion in the Stark problem is bound whenever E < and 
\S\ — S ^ S crit (E) = E 2 /(GM), and it is nearly-Keplerian when S <C S crit - In this latter 
case the orbit consists of a precessing ellipse of varying eccentricity and inclination. In the 
2-D case [z = p z = S z = 0, say) the maximum eccentricity reached, e max , is always unity, 
and thus these orbits momentarily become radial (at which time the circulation of the orbit 
changes from prograde to retrograde, or vice versa); the minimum eccentricity, e m ; n , normally 
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occurs where the line of apsides is parallel to the Stark vector and can lie anywhere between 
and 1. In the 3-D case, conservation of the angular momentum component forces e max < 1, 
and now the inclination also varies between minimum and maximum values. In all cases, 
e m i n and e max are constants of the motion; this conveniently allows fine control over the range 
of eccentricities encountered during testing, regardless of the length of the integration. An 
example of nearly-Keplerian (S = 0.12S , cr i t ), 2-D motion is given in Figure |l|. 



2.2. Behavior of the Wisdom-Holman Mapping 

Using the WH mapping as a numerical integrator for the Stark problem is equivalent 
(within round-off error) to replacing the Hamiltonian ([[J) with a "nearby" mapping Hamil- 
tonian, and solving the resulting equations of motion exactly (see [Wisdom fc Holman 1991| ). 
In the present case, a second order mapping Hamiltonian corresponding to equation ([!]) is 

(v 2 GM\ 

H^=l*L- — \- 2n5 27T (Qt - tt)(S • x), (2) 

where ^(x) is a periodic delta function (with period 2tt), Q = 2k / At is the mapping 
frequency, and At is the associated integration step size. (One physical realization of this 
Hamiltonian is an orbiting satellite performing periodic, short-duration burns of its engine.) 
The mapping Hamiltonian differs from the original by a series of high frequency (Q and higher 
harmonic) terms. In general, the long-term evolution under H m&p can be expected to remain 
"close" to the true evolution as long as the mapping frequency exceeds all fundamental 
dynamical frequencies in the problem; this is known as the averaging principle. The second 
order integration step corresponding to H mSuP is 

I{x, p; At, S) = D(x, p, At/2) K(S, At) D(x, p, At/2), (3) 

where D represents a drift along an unperturbed Keplerian orbit and K represents a mo- 
mentum kick due to the perturbation S (x is left unchanged by the map K). We remark 
that the Kepler step D is most efficiently computed using the Gauss / and g functions (e.g., 
Danby 1992| ). 

A typical example of the evolution under if map is shown in Figure ||, which plots the 
fractional energy error committed for a pair of 2-D integrations using 100 (dashed curve) 
and 1000 (solid curve) points per orbit; the beginning orbit, similar to Fig. [1|, had an initial 
eccentricity eo = 0.9 and a Stark perturbation S = 4 x 10 _3 S cr it directed at a 45° angle from 
the initial line of apsides. In contrast to the bounded, oscillatory energy errors typically 
observed in symplectic integrations, in this case the errors undergo a random walk towards 
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Fig. 1. — An example of bound, nearly-Keplerian Stark motion; the Stark vector is directed 
along the +x direction. The initially prograde, moderate eccentricity (eo = 0.8) orbit begins 
with its line of apsides along the x-axis. The orbit subsequently suffers a clockwise precession 
and becomes increasingly radial. Upon reaching e = 1, the orbit switches to retrograde 
motion and counterclockwise precession (not shown); the orbit eventually returns to prograde 
circulation in the lower left quadrant, and the cycle continues. 
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unity, even using 1000 points per orbit. Eventually, in fact, these numerical orbits become 
unbound and escape to infinity — a qualitatively incorrect result for the long-term motion of 
the particle. This is suggestive (though not conclusive) evidence for numerical chaos, which 
if present must arise from the integrator itself since the analytic motion is fully integrable. 
Although this result does not, strictly speaking, violate the averaging principle (because 
the mapping frequency does not exceed the pericentric passage frequency when the orbit 
is nearly radial), it is nonetheless surprising that the orbits are so unstable under such a 
modest perturbation — particularly since the relative perturbation strength at periapse is 
even smaller. 

A less typical (but not infrequent) example of the evolution under H map is shown in 
Figure |3], which is analogous to Figure |2] except that the Stark vector S = 4 x lO -5 ^^ 
and is parallel to the initial semi-major axis. As before, the 1000-point-per-orbit integration 
initially shows bounded, oscillatory energy error which degrades significantly once periapse 
in no longer resolved. In this case, however, both integrations appear to possess bounded 
energy errors at late times, albeit at a much higher level than is exhibited initially. This 
long-term stability generally persists indefinitely; we have let particular simulations run for 
~ 10 9 orbital periods with no apparent change in the error bound. 

A strongly suggestive illustration of the underlying distinction between the two types 
of behavior is show in Figure |4]. The figure contains three panels, each showing the regions 
of stability in the At — S plane for identical initial conditions aside from a change in orbital 
phase; in each case eo = 0.9 and the Stark vector is parallel to the initial line of apsides. 
The initial mean anomalies of the orbits, from top to bottom, are Aio = 0, 2n/3, it. The 
'islands' of stability present in the panels display a clear pattern: they are all centered on 
step sizes which are a rational fraction of the orbital period, and are most prominent at step 
sizes corresponding to an integral number of points per orbit. This is clear evidence that 
the stable behavior is the result of step size resonances, and that otherwise the mapping is 
generically unstable; a formal analytic analysis supporting this assertion will be presented in 
the following section. Note also how the size of any individual island varies with the initial 
orbital phase — its area reaches a maximum when one of the integration steps regularly lands 
near periapse. The 2:1 (At/t or b = 1/2) resonance island, for instance is large in the Aio = 
and A4o = Ti panels (for which every other step — the odd or even ones, respectively — falls 
near pericenter), but small in the M.q = 2ir/3 panel (where all steps fall rather far from 
pericenter). This points to another necessary condition for stability; namely, stability does 
not require that periapse be resolved (which is not even possible in two dimensions, since 
e ma x = 1), but merely that it be sampled. We explain these observations in detail in § [2~3 . 



Although true dynamical chaos is the most natural explanation for the unstable be- 
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log t/t orb 

Fig. 2. — The fractional energy error over time committed in the integration of an S = 
4 x 10" 3 S' cr it, eo = 0.9, 2-D Stark orbit using the Wisdom-Holman method; the solid curve 
is for At = lCT 3 t or b, the dashed for At = 10 _2 t orb . The initially oscillatory behavior of the 
solid curve fails when e increases to where the step size no longer resolves periapse; at late 
times, both orbits random walk in energy and eventually go unbound. Note that although 
the initial separation of the curves is just what would be expected from a second order 
integrator, the separation is not even qualitatively maintained over the long run. 
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log t/t, 



Fig. 3. — As Figure |2|, but for S = 4 x 10 5 S cr i t (with S parallel to the initial line of apsides). 
Note the bounded energy error at late times. 
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Fig. 4. — Islands of stability in the At — S plane (regions smaller than a few pixels are 
noise). Each panel corresponds to identical initial conditions except for the starting orbital 
phase. The initial mean anomalies, from top to bottom, are Aio = 0, 2-7r/3,7r. The islands 
all center on rational fractions of At/t OT b, but their shapes depend strongly on Mo. See 
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Fig. 5. — The Lyapunov times £l across a horizontal slice (at S/S cr i t = 0.04) of the Ai Q = n 
panel of Figure §; values of ^ 10 5 t O rb should be regarded as lower limits due to the finite 
length of the integration. The location of the stable islands is clearly visible, indicating that 
these are regions of true dynamical regularity (and not bounded chaos). 
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havior present in Figure 0, none of the evidence presented so far proves that the ostensibly 
stable trajectories in Figure |3| are in fact dynamically stable — i.e., regular. To test this, 
we computed the Lyapunov times (e.g., [Bennetin, Galgani, fc Strelcyn 197"S| ) for a series of 



trajectories corresponding to a horizontal slice (at 5* = 0.04S , cr it) across the bottom panel 
of Figure [§ The result is shown in Figure [|; note that all values of £l ^ 10 5 t or b should be 
regarded as lower limits due to the finite length of the integration. The figure provides con- 
vincing evidence that the 'stable' islands are, in fact, dynamically stable, and not an example 
of bounded chaos; off the islands, by contrast, the motion is strongly chaotic. Examination 
of sample surfaces of section confirmed this result, and also verified that the integrated mo- 
tion inside the stable islands remains close (in terms of explored phase space volume) to the 
analytic trajectory for arbitrarily long times. 

The preceding results do not qualitatively change when the general three-dimensional 
problem is considered. The stable islands are still present, although they shrink — and even- 
tually disappear (except the one for which At resolves e max ) — as the Stark vector is rotated 
from the x-axis (the initial line of apsides) to the z-axis (perpendicular to the initial orbital 
plane); however, the same phenomenon also occurs in the 2-D case when S is rotated from 
the x-axis to the y-axis. In both cases the disappearance of the islands coincides with the 
condition 

Cmin ~ 0; the connection is explained in § |2.3.2| . 

Finally, a few words about numerics. Because of the extreme range of orbital eccen- 
tricities encountered, the Stark problem places severe stress on the Kepler stepper (the drift 
operator in equation ||) used in the integration; the kick operator, by contrast, is completely 
trivial. In fact, some of the stepper routines at our disposal (originally written with low to 
moderate eccentricities in mind) failed outright for nearly radial orbits, particularly at large 
step sizes. Among those proving robust, however — including one using extended precision 
arithmetic throughout — the choice of stepper did not alter the stability of the mapping aside 
from minor changes in the precise boundaries between the stable and unstable regions. Our 
testing indicates that straightforward implementation of the universal variables formulation 
of the Kepler stepper, combined with constrained solution of Kepler's Equation (i.e., one 
guaranteed to converge) and renormalization of the final radius and/or velocity vector (to 
explicitly enforce energy conservation), produces a nearly optimal routine in terms of both 
efficiency and robustness. 



2.3. Non-Linear Stability Analysis 



Analytic examination of the non-linear stability of the WH method was carried out using 
the resonance overlap criterion described in Wisdom & Holman (1992). In brief, the method 
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involves expansion of the mapping Hamiltonian as a sum of resonant terms, with only the 
term representing the step size resonance under consideration being retained. Expansion 
of this Hamiltonian about the resonant value of the canonical momentum then leads to 
expressions for the width and libration frequency of that particular resonance. Overlap, 
and subsequent instability, occurs when the allowable libration amplitudes (in the energy 
oscillations, say) of adjacent resonances exceeds their separation. 



2.3.1. Hamiltonian Development 

To begin the analysis we rewrite the mapping Hamiltonian in equation (0) to explic- 
itly show the interaction of the step size dependent terms with the terms of the original 
Hamiltonian. Using the Fourier representation of the delta functions, 

J oo 

5 2 n(nt) = — E cos(ifit), (4) 

2,71 ■ — 

l= — CO 

the first order mapping Hamiltonian can be written as 

/ GM \ 00 
Hmap = [p 2 J - Yl cos(iQt)S ■ X. (5) 

^ % — — CO 

The time step, or interval between delta functions, is At = 2n/Q. To simplify the analysis 
we here consider the 2-D case with the Stark vector parallel to the x-axis. The quantity S ■ x 
is then S x x = S x rcos(9) = S x rcos(ui + /). Here r is the "heliocentric" distance, u>i is the 
longitude of pericenter, and / is the true anomaly. Thus, the Hamiltonian is 



CO 



Hmap = (p 2 — -^-] ~ coB(iQt)S x r cos(ui + f) (6) 



%— — CO 



GM\ 00 

p 2 J — E cos(i£lt)S x (r cos/cosu;/ — r sin/sinu;/) . (7) 



Next, we expand r cos / and r sin / in terms of Bessel functions, assuming a Keplerian orbit. 

r 3e 00 1 

-cos/ = cos(£)-e = -— + 2Y,TJ'k(ke)cos{kM), (9) 
a 2 fe=1 k 

r , , 00 1 

-sin/ = Vl - e 2 sin(E) = y/l - e 2 ^ -rJk-i(ke) sm(kM), (10) 
a k=1 k 

where a is the semi-major axis, E is the eccentric anomaly, e is the eccentricity, Ai is the 
mean anomaly, k is an integer sequence, and J^-iike) and J' k {ke) are a Bessel functions and 
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their derivatives flBrouwer fc Clemence 196T ). Thus, the Hamiltonian becomes 

GM\ 



H, 



map 



P~ 



— dS x cos X cos(zf2t) 

i=— oo 

+aS x sin toi ^ cos(if2t) 



00 1 

-f + 2 E^(Mcos(A;>l) 
fc=i 

00 

Vl - e 2 -rJk-i(ke) sin(kM) 
k=i k 



(11) 
(12) 

(13) 



where we ignore for the moment that e and A4 are not appropriate canonical variables. The 
trigonometric factors can be re-arranged and combined: 



map 



GM\ 



1 



(14) 



+aS x cos u>i cos(iQt) 



3e 



%=— 00 
00 00 



— aS x cosui E E —J' k {ke) cos(/c.M — iVLt) 
i=o k=i k 

00 00 

+aS x sin u)\ J2 E v / T^^7 Jfc-i(fce) sin(fc.M - iOt). 

i=0 k=l k 

Next, we make two related simplifying assumptions. The first is that S x is sufficiently small 
that the timescale for changes in the eccentricity e and the longitude of pericenter lo\ is 
much longer than the orbital period or other timescales in the system. Thus, e and u)\ 
will be considered as constants. By this assumption the terms in the first summation are 
strictly time dependent and will not contribute to the equations of motion resulting from the 
Hamiltonian. We thus ignore those terms. The second assumption is that the eccentricity is 
sufficiently large that the \/\ — e 2 factor is negligibly small. Dropping the final summation, 



H 



map 



P~ 



GM\ 



-) 



(15) 



1 



k 

i=0 k=l " 

Re-writing the Hamiltonian in canonical variables: 



— aS x cosui E X/ T^fc(^ e ) c °s(fc7W — iQt). 



H 



map 



(GMf 
2L 2 



aS x cosui E X/ T^fc(^ e ) cos [^(^ — Ui) — iQt], 
i=o k=i k 



(16) 



where the momentum L = v GMa is canonically conjugate to the mean longitude A = 
uj 1 + Ai. The arguments to the cosine terms are now clearly the "step size resonances" in 
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the Hamiltonian. Such a resonance occurs when one of the terms of the form k(X — u>i) — iVLt 
is slowly varying. Given our assumption that uj\ is roughly constant, this happens when 
kX — ifl PS 0, i.e., when the orbital period is rationally related to the step size. 

Next we make a canonical change of variables that focuses on one of these terms. For 
this we use the mixed-variable generating function, 



F 2 = (X-ui- -Qt)E + AA, 
k 



which results in the following transformation: 

dF 2 



a 
X' 
L 



<9E 
dF 2 

dA 
dF 2 

dX 



= (A - LOl) - -fit, 
= A, 

= A + E. 



(17) 



(18) 



The new Hamiltonian is 
H' = 



H + 



dF 2 



dt 
{GMf 



2(A + S) 2 

where only the resonant term has been retained. 



1 % 
— aS x coscj/— J' k (ke) cos(ka) — — fi£, 



(19) 
(20) 



Since we expect the canonical momentum to be constrained to a narrow range of values 
at resonance, E should vary little from the resonant value. Consider only the momentum 
terms in the Hamiltonian: 



H'o 



(gmy 



ny 

2(A + E) 2 k 

which can be expanded about the resonant value of E (S*): 

1 d 2 H' 



-n - -"ols* + 



as 



;e-e*) + 



2 as 2 



*\2 



^E - E*) 



(21) 



(22) 



The first term is a constant and can be ignored. The second term defines the resonance 
condition 



dH> 



as 



(GMf 



- -n 



n* - -n 

k 



o, 



(23) 



i: (A + E*) 3 k 

which is identical to the earlier statement of what constitutes a step size resonance (where 
n* is the corresponding mean motion). This leaves only the quadratic term, 

{GMf 

(A + £*) 4 



7 



d 2 K 



as 2 



-3 



-3- 



1 



(24) 
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The Hamiltonian can now be reduced to a standard form: 




(25) 



where AS = E — E* and (3 = — a*S x ( J' k (ke)/k) cosc^. We have also assumed that the 
variations in semimajor axis are small enough that we are justified in using the resonant 
semimajor axis, a*, in f3. This is just a modified pendulum Hamiltonian. For (3 < 0, the 
argument of the cosine, ka, oscillates about the values 2irm (where m is an integer from 
to k - 1). 

The half-width of the resonance is normally defined by the largest value of the momen- 
tum for which the trajectory still librates rather than circulates. This is 



It is important to note that the resonance width depends upon k but not upon i. Thus, in 
the vicinity of a particular step size the closest k : 1 resonance will be most important. From 
the definition of E, AL = AS. And for small amplitude oscillations, Aa/a ~ 2AL/L and 
An/n ~ —3AL/L. Following the pendulum approximation further, the frequency of small 
oscillations (the libration frequency) is given by 



In the following section we use the components of the Hamiltonian development to support 
a geometric interpretation of the resonance overlap criterion. 



The libration about resonance derived above manifests itself geometrically in terms 
of an oscillation about pericenter of the orbital mean anomaly associated with successive 
integration points. As an aid in visualization, consider, for example, a step size equal to 
the mean orbital period. If the integration is started at periapse, then each subsequent 
step will also begin (and end) near periapse, aside from a small drift induced by the Stark 
perturbation. In stable libration, however, this drift away from periapse is replaced by 
an oscillation centered on periapse; in this case, the particle is stably trapped in the 1:1 
resonance. Figure |6| shows an example of libration in the 3:1 resonance, in which each of the 
3 points per orbit oscillates around a fixed value of the mean anomaly M. (namely M. = 
—2tt/3, 0, 2n/3). When this oscillation is stable, as it is in the figure, the behavior shown 




(26) 




a* 



(27) 



2.3.2. Geometric Interpretation 
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Fig. 6. — An example of stable libration of the mean anomaly in a 3:1 resonance. The curves 
display two periods, the shorter, i^, being the primary libration period and the longer, t ccc , 
being the timescale for the eccentricity to vary between e m i n and e max . Libration is stable 
only when tub ^ t ccc and its total amplitude is ^ 2n/k (for a k : i resonance). 
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in Figure |3] results; otherwise, the behavior is that of Figure Q The reason is that the stable 
oscillation systematically cancels out the energy errors that would randomly accumulate in 
its absence, thereby stabilizing the long-term motion (cf. the dashed curve in Figure ^ 
where the initially linear error growth turns into an oscillation on a timescale corresponding 
to the libration period). Also note that this libration cannot remain stable as e — > 0, since 
pericenter becomes undefined there; this qualitatively explains the disappearance of the 
stable islands for e m ; n ~ which was mentioned in § |2.2j . More formally, since the width of 
the k : 1 resonance AS oc [J' k {ke)] l l 2 (see equation |^), it follows immediately that AS — > 
as e — > 0. 

With this simple picture in mind a scaling argument describing the shape of the stable 
islands (Figure |]) is easily obtained. To do this, we first define the resonant timestep 
corresponding to the k : 1 resonance, At res = t OI h/k. For a nearby integration with timestep 
At ~ At res , define next a "drift" timescale, tdriftj representing the time needed for the mean 
anomaly of every k th integration point to drift by 2n/k (due to the small mismatch between 
At and At res ): 1/idrift = |l/At res — 1/At\. In qualitative terms the libration will go unstable 
when t drift ~ t lib , where ti ih ~ (S/ Scnt)^ 2 t olh can be derived from the results of the previous 
section. This can be rewritten 

At 2 ( S \ 1/2 

| At - At res | ~ — - -— , (28) 



^orb V S, 

where \At — At res | <C At res and S <C S crit are implied. Since At ps to^/k, it follows that 
the width of the islands should roughly scale like k~ 2 ; in addition, their boundaries should 
approximate a parabola in the neighborhood of At res (S oc (At — At res ) 2 there). To the 
extent that they can be reliably measured, these scalings hold quantitatively in Figure ^ 
(the bottom panel is cleanest in this regard). A surprising corollary of this is that decreasing 
the perturbation at fixed step size can destabilize the calculation by taking it outside the 
local stable island. 

We can also use the results of the previous section to estimate S mSjX (k), the maxi- 
mum height of the k : 1 resonance island. In this case we use the relations AL/L ~ 
[(J{,(fce)/fc)(S'/S' cr it)] 1//2 and Aa/a ~ 2AL/L (assuming small amplitude oscillations) and 
note that the libration in a, the semi-major axis, cannot remain stable when the k : 1 
resonance of this orbit overlaps the k + 1 : 1 resonance of the nearby orbit with semi- 
major axis a + Aa. The latter condition occurs when Atk(a) ~ Atk+i(a + Aa) (where 
Aifc(a) = t OI b(a)/k), which implies Aa/a ~ k^ 1 for S ~ S^x. Equating the two then gives 

s - {k) ~ [mm) «* < 29 > 
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To order of magnitude we can take J' k (ke) ~ [(1 - e^/V^k] exp (-(1 - e 2 f' 2 k/i) (e.g., 
Abramowitz fc Stegun 1968| ; the formula becomes exact for e ~ 1 as k — > oo). Thus 
when k <C k CTit = 3(1 — e 2 )~ 3//2 (i.e., when At does not resolve periapse), S max oc A; -1 / 2 ; 
this is close to the actual scaling in the bottom panel of Figure ^ (for which the empirical 
values of S max are least ambiguous). Note that the estimated S max reaches a minimum for 
k ~ fccrit and increases exponentially for k ^ k crit . This suggests that the mapping should be 
stable whenever periapse is well resolved; all the numerical testing we have done confirms 
this prediction. In Figure f| this phenomenon does not occur since e max = 1 and hence 
no fixed step size can resolve every periapse; in the 3-D case, where e max < 1, a 'wall' of 
stability is created enclosing the entire region At < t OT h/k cr it(e mauX ). Even though the relative 
perturbation strength is often minuscule at periapse, it therefore appears the stability of the 
WH mapping requires that it be well resolved regardless. The concomitant loss of efficiency 
for highly eccentric orbits is obviously enormous. This motivates the search for more robust 
methods not subject to this limitation, the topic to which we now turn. 



3. Modified Wisdom-Holman Mappings 

There are at least two important advantages to integration methods whose stability 
hinges only on resolution of the highest frequencies associated with the perturbation forces, 
expert, instead of those intrinsic to the unperturbed motion (which we presume the method 
to handle exactly), uj ot \,. The first, of course, is efficiency: if cj rb ^> expert, such a method 
can use a much larger timestep — up to uj orh /uj pen times larger — than a scheme that must 
explicitly resolve uj ot ^. A more subtle advantage concerns the unavoidable loss of energy 
accuracy due to round-off errors near periapse. More specifically, if the motion is nearly- 
Keplerian with eccentricity e no algorithm based on Cartesian phase space coordinates, that 
also samples pericenter, can maintain better than N + log |1 — e| digits of energy accuracy, 
where N is the number of arithmetic digits carried. The proof follows immediately from 
E = p 2 /2 — GM/r = (rp 2 — 2GM)/{2r) and the fact that rp 2 « (1 + e)GM near pericenter. 
Note that doing selected intermediate calculations in extended precision arithmetic will not 
improve on this; the mere act of rounding the output values of x and p to N digits is 
sufficient to do the damage. In principle this problem can be circumvented by recasting the 
dynamics in terms of osculating orbital elements (for example), but the formulation is likely 
to be awkward (and probably inefficient) since the equations are not being orbit-averaged. 

The search for WH-like algorithms that are robust in the above sense is thus well moti- 
vated. In the following sections we will limit ourselves to examining two recently proposed 
variants of the basic WH scheme, both still exact for unperturbed Keplerian motion. The 
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utility of a more divergent line of methods, based on Stark instead of Kepler motion, will be 
discussed in § (|. 



3.1. The Regularized Wisdom-Holman Mapping 

Regularization of the WH mapping through the use of an extended phase space has been 
proposed by Mikkola (1997). In this approach the integration substitutes a regularized time 
s for the physical time t, the defining relation between the two (as in K-S regularization; 
e.g., [Stiefel fc Scheifele 1971[ ) being ds = dt/r; further replacement of x and p by the 



corresponding K-S variables is straightforward but entirely optional, and is not considered 
here. The constant steps in s used by the method naturally sample pericenter more densely 
than constant t-steps would (although it is still not quite resolved — that would require ds ~ 
dt/r 3 ^ 2 ), offering hope of increased reliability at high eccentricities. The substitution of s for 
t is done by extending the phase space of the original Hamiltonian to include t and E (the 
total energy) as an additional pair of conjugate coordinates; for details, see Mikkola (1997). 

The numerical performance of this method when applied to the Stark problem is illus- 
trated in Figure |^, which plots the errors in energy, angular momentum, and position of the 
method (using 100 points per orbit) relative to the analytic solution for parameters identical 
to those in Figure |2|. Accuracy and stability in this case are excellent, as they were in every 
case we tried (see § |5| for additional examples). Examination of the corresponding surfaces 
of section confirmed that the integrated motion was regular, even for integrations using 
fewer than 10 points per orbit — a remarkable feat considering the eccentricities involved! It 
is in fact possible that the corresponding mapping Hamiltonian is integrable for arbitrary 
timesteps, although we have not proven this. 

This result demonstrates not only that methods which are stable for timesteps At ~ 
1/^pert ^ 1/^orb exist, but also that they can be simple and efficient. In fact the regularized 
WH map for this problem is noticeably faster, per step, than the original mapping (cf. 
Table p] in § ||); this is because the regularized Kepler stepper does not need to solve Kepler's 
Equation, and hence is much faster than the original one (this applies only to problems of 
the perturbed two-body type). Because of this, it could even be argued that the regularized 
method is always superior to the original for this class of problems, even when eccentricities 
are low — although the speed increase will be noticeable only when the cost of the Kepler step 
is a significant fraction of the total. On its own, however, this method cannot handle close 
encounters with perturbing objects; a promising strategy for incorporating this capability 
into either the original or the regularized method is the subject of the following section. 
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Fig. 7. — A typical example of the performance of the regularized Wisdom-Holman mapping 
(§ |3.ip applied to the Stark problem. The relative errors in energy (solid), angular momentum 
(dotted), and position (dashed) are plotted for an integration utilizing 100 points per orbit 
and the same initial conditions and Stark vector as Figure Unlike the original mapping, 
the modified method is completely stable — the energy error is bounded and the lag in orbital 
phase grows linearly with time. 



-21 - 



3.2. Potential-Splitting Methods 

A well-known limitation of the WH method (and related symplectic schemes, such 
as leap frog) is that the timestep cannot easily be varied without destroying its desirable 



symplectic properties, such as long-term energy conservation (e.g., |Gladman, Duncan, fc 



|Candy 1991| ; |Skeel fc Gear 1992]) . This makes the creation of adaptive symplectic algorithms 



capable of handling close encounters a delicate undertaking. Duncan, Levison, & Lee (1997; 
see also Lee et al. 1997), bulding upon the approach of Skeel & Biesiadecki (1994), have 
recently developed a symplectic, multiple-timestep generalization of the original WH method 
which adds the ability to resolve close encounters without seriously compromising its overall 
efficiency We will refer to it as the 'potential-splitting' (PS) method because of the way it 
splits the potential of each perturber into a series of radial zones centered around it. We 
find the approach interesting not only for its versatility in handling close encounters within 
a symplectic framework, but also because the scheme is amenable to regularization. To our 
knowledge this latter possibility has not been explored elsewhere; our subsequent discussion 
of the PS approach will concentrate on determining the utility of such a merger. 

In brief, the PS method works as follows. Consider for simplicity a two-body orbit 
perturbed by a single point mass m at a fixed position x p , which generates a potential 
U(x) = —Gm/\x — ajp|; the dominant central mass, M, is assumed to be at the origin. 
This is nothing but the two fixed point problem that will be used in § |5.1| for comparative 
testing. The PS method divides the potential U — or more conveniently, the perturbation 
force F = —VU — into a series of shells, F = YljfLo Fji eacn ( exce Pt ^o) being non-zero 
over only a finite range of p — \x — x p \. Introducing an ordered sequence of radii pi (i > — 1), 
where p_i = oo and Pi/pi-i = const. < 1 (z > 0), one then defines Fj = hj(y)F, where 
y(x) = (p(x) - pi)/(pi_i - pi) (subject to < y < 1, implying p { < p(x) < p;_i) and 

Pj < P< Pj-i, 

h j{y) = { Pj+i < p< pf, (so) 

otherwise, 

where the splitting kernel n(y) is a monotonic function satisfying k(0) = and k(1) = 1. 
(The preceding are not the most general definitions, but are the ones we will use.) As 
noted by Lee et al. (1997), it is also desirable for the derivatives of k to vanish at the 
endpoints as it increases the smoothness of the transition between neighboring Fj] they 
suggest n(y) = y 2 (3 — 2y), which is the unique cubic satisfying k(0) = 0, /t(l) = 1, and 
k'{0) = k'{1) = 0. 

The algorithm proceeds by splitting the base timestep At into an integral number of 
smaller steps whenever the test particle enters a zone interior to the one it previously resided 
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in — i.e., whenever it happens to approach the perturber sufficiently more closely. At each 
such subdivision, a specific kick component Fj is applied in such a way as to allow further 
subdivisions if necessary yet keep the entire process symplectic; eventually all higher num- 
bered Fj vanish, the recursion is terminated, and the intervening Kepler drift about M is 
performed. The net effect is to take many small steps during close encounters; far from the 
perturber, the method reduces to the standard WH mapping (albeit in kick-drift-kick form 
instead of drift-kick-drift form; cf. equation ||). For further details, consult the references. 

We now propose several enhancements to the basic PS method which we have found 
can significantly improve its stability. The first is a modified kernel function, n(y). The 
motivation is simply to produce a fully analytic decomposition of F, meaning one for which 
all derivatives of n{y) vanish at the endpoints, instead of merely the first; this is conceptually 
similar to the creation of a bump function." Many forms are possible; the one we have 
found most useful is 

(31) 

Although more costly to evaluate than the original polynomial (the new code was about 
30% slower), the resulting method was stable for larger At and longer periods of time than 
the original in every test tried (including the two fixed center problem, § |5.1|) . On the other 
hand, since in many cases the polynomial kernel performed nearly as well, we do not claim 
our modified kernel is uniformly superior; we have, however, found it a useful alternative in 
situations where the original seems to behave poorly. 

The second enhancement involves force-center switching. Whereas the basic PS method 
takes Kepler drifts about M regardless of the test particle's proximity to the perturber m, 
it is clearly advantageous to execute drifts about m during very close encounters. This 
capability can be neatly incorporated into the PS scheme as follows. Note first that the 
potential U is not split at all when p > p ; this is the regime in which F = F , where no 
subdividing is done and the method is equivalent to the usual WH approach. Now define a 
value j max such that F = Fj m ^ whenever p < Pj max ; this implies that all Fj are identically 
zero for j > j max - Thus when p < Pj max the full Hamiltonian is again available and can be 
separated into a piece representing an m-drift and a piece representing the "perturbation" M 
with potential V = — GM/\x\. Next, modify the rules for subdivision so that when p < Pj max , 
kick(l//2)-drift(m)-kick(y/2) is done instead of Mck(LT, mM /2)-drift(M)-kick(C/ imaj£ /2); since 
Uj max = U here both forms are equally valid. The resulting method now uses drifts around 
m (without further subdivision of the timestep) whenever the encounter is close enough. 
Although any j max > may be used, it should obviously be chosen so that m dominates 
the dynamics for p < Pj max . It may appear that switching splittings in this fashion breaks 
the symplecticity, and this is entirely possible; if true, however, we have not found it to be 



K(y) = — < 1 + tanh 



2y-l 
y(i-y) 
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a problem. In tests with the two fixed center problem, for example, energy conservation 
remained stable even after thousands of switchings. This may be due partly to the fact that 
the local timestep is very small when the switching takes place, minimizing any systematic 
errors it may commit. The technique appears a promising one in any event. 

The final enhancement we considered was the incorporation of regularization into the PS 
framework. This is in fact quite straightforward to accomplish; one simply splits the regular- 
ized Hamiltonian instead of the original one — the Fj becoming regularized force components, 
and so on. The usefulness of doing this is also clear: whereas the basic PS method shares the 
instability of the WH mapping it is based on when orbits are eccentric, the regularized PS 
method remains robust here. We have verified this directly in the case of the Stark problem; 
in particular, the instability in the unmodified PS method was found to persist even when 
splitting of the Stark potential at small radii was included — not to mention the fact that 
doing this made the method extremely inefficient! As we will clearly demonstrate in § |5], the 
regularized PS algorithm appears completely robust in this regard. 



4. Stark-based Integration Schemes 

A more radical strategy for creating integrators reliable in the way outlined in § § is 
not to make minor transformations to the original Hamiltonian splitting, but instead to 
rethink the Kepler splitting entirely. In this case the primary motivation for altering the 
splitting is not to produce a faster method, as was true for the WH mapping (compared 
with leap frog, say), but rather one that is more stable. The challenge is that each piece 
of the new Hamiltonian splitting should be integrable and efficiently soluble if the mapping 
is to be practical. Although integrable problems in general are a precious commodity, the 
preceding analysis provides two obvious candidates: the Stark problem and the two fixed 
point problem. (The two are in fact closely related, since the latter reduces to the Stark 
problem in the limit |aj p | oc m 1 / 2 — > oo, where m and x p are the mass and position of 
one of the fixed points.) In this paper, however, we will only consider mappings based on 
a Stark splitting of the Hamiltonian, arguably the simpler of the two (it having one less 
free parameter). For comments on the use of the two fixed point problem as the basis of a 
symplectic mapping, refer to § |^. 



4.1. The Time- Reversible Stark Method 



The first Stark-based method considered was derived from the original WH mapping 
by simply replacing the Kepler stepper in that method with a Stark stepper, where the 
value of the Stark vector for a given step was self-consistent ly taken equal to the value of 
the perturbation force at midstep; consequentially, the intervening 'kick' vanishes. (Hence 
computation of the Stark vector requires iteration, at least in principle.) An integration 
therefore consists solely of a sequence of Stark steps, the Stark vector for any particular 
step representing a local best-fit to the perturbation. Note that it would be misleading to 
think of this method as symplectic. The basic problem is that in each step the Stark vector 
depends on the current position of the test particle, yet this coordinate dependence in the 
mapping Hamiltonian is not accounted for by the Stark step itself. Thus although it is 
perfectly reasonable to think of any one step as being symplectic — since a strictly constant 
Stark vector which cancels the corresponding kick for that step is easily manufactured — the 
unconstrained manner in which the vector changes destroys the self-consistency required 
for the sequence of successive steps to exhibit symplectic behavior. The situation is closely 
analogous to attempting to vary the step size in, for example, leap frog or the WH method: 
although any one step is obviously symplectic, successive steps are not part of the same 
Hamiltonian flow and hence the integration does not display coherent symplectic properties 
such as long-term energy conservation. By construction, however, the method is explicitly 
time-symmetric and hence might still exhibit good energy behavior; we will refer to it as 
the time-reversible Stark method. Also note that for eccentric orbits the perturbations near 
pericenter will generally be nearly stationary, so that the Stark approximation should be 
an excellent one here — giving hope that the scheme will maintain stability in such cases. 
Disappointingly, this hope was quickly dashed by the results of our numerical tests; for 
details, see § |. 



4.2. Regularized Stark Mappings 

The regularization and extended phase space techniques employed by Mikkola (1997) can 
also be harnessed to create a regularized, fully symplectic Stark-based mapping. Consider in 
particular a time-dependent, perturbed two-body Hamiltonian of the form H = p 2 /2 — l/r + 
U(x,t) (where r = \x\), and introduce a fictitious Stark potential —S(t) ■ x that reproduces 
U(x,t) as closely as possible; the critical restriction is that S can depend only on t, not on 
x. For example, if U(x,t) = —Gm/\x — x p (t)\, then S(t) = Gmx p (t)/\x p (t)\ 3 would be 
optimal whenever |je| <C \x p (t)\ (and otherwise the perturbation is not Stark-like at all, so 
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using a Stark splitting would be pointless to begin with!) Now rewrite H as 

H(t) = ^----S(t)-x + 5U(x,t), (32) 
2 r 

where SU(x,t) = U(x,t) + S(t) ■ x. Letting £ = —H(0) be the total binding energy and 



extending phase space using ds = dt/r (cf. § |3.1| ), the extended, regularized Hamiltonian is 

+ r[£ + 8U(x,t)\. (33) 



H 



p 2 1 

2 r 



The grouping is intended to show how H = H + Hi is to be split. The first piece, H , is 
clearly just the (regularized) Hamiltonian for the time-dependent Stark problem; however, 
since Ho is independent of £ Hamilton's equations imply that the coordinate t is a constant 
here (it is advanced only by the second piece of the Hamiltonian, Hi). Therefore, S(t) 
is a constant vector for the duration of the H step, and we have succeeded in creating 
an obviously symplectic algorithm based on perturbed Stark motion instead of perturbed 
Kepler motion. In addition, the (optional) regularization used can be expected to provide 
the same protection against instability as occurred in the regularized WH method. 

Although we came upon it independently, the above method (without regularization) 
is very similar to the one created by Newman et al. (1997; see also |Grazier 1997|) . Their 
motivation, by contrast, was in exploring its use in treating close encounters; their approach 
(as we understand it) represents a modification of the method of Levison & Duncan (1994) 
in which the Stark splitting is employed only during close encounters, just after switching 
force centers to the nearby perturber. Irrespective of motivation, it must be kept in mind 
that there is a significant price to pay for using a Stark splitting, since Stark steps (in our 
experience) are roughly 30 times slower than the Kepler steps they are replacing. Unless 
the force calculation strongly dominates the execution time, this implies that an increase 
in the step size by a factor of 30 must be permissible for the change to pay off. This is 
clearly a severe constraint! On the other hand, if the method is stable for At ~ l/u pert , 
and uj orh /u pcrt ^ 30 (cf. § |3|), then such huge improvements would be plausible. When the 
original Kepler mapping is also stable, however, it is quite unlikely gains of that magnitude 
would be possible — particularly since the use of symplectic correctors ( |Wisdom, Holman 



fc Touma 1996| ), which we have not taken advantage of, can often substantially boost the 



accuracy of the original mapping with no loss in overall efficiency 



5. Comparative Simulations 
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5.1. The Two Fixed Point Problem 

Having described four alternatives to the original WH scheme — the regularized WH 
mapping, the PS method (with and without our purported enhancements), and two types 
of Stark-based schemes — we now wish to provide some practical insight into their relative 
performance. In this section the two fixed point problem, whose dynamics is understood in 
complete detail, is used as a test problem; in § |5.2| a more generic test problem from the area 
of galactic dynamics is used to estimate performance under more "typical" conditions. 



5.1.1. Orbital Motion 



The two fixed point problem (e.g., Pars 1965 ) represents the motion of a test particle 



in the field of two gravitating point masses, mi and rri2, held at fixed positions, X\ and x 2 , 
the Hamiltonian per unit test mass is 

2 | a; — x\\ \x — x 2 \ 

Like the Stark problem, it is fully integrable and possesses three constants of motion; as noted 
earlier, it in fact reduces to the Stark problem in the limit (for example) |a?2| oc m 2 — ► oo. 
The problem is separable in confocal coordinates and can be solved analytically in terms of 
elliptic functions and integrals. The three constants of motion are the energy E, the angular 
momentum component along the direction x 2 — Xi, and a separation constant a (an explicit 
formula for which can be found in Lessnick (1996)). 

Although the problem is rather artificial (we can think of no good physical analogies), 
it does allow close encounters to be introduced in a controlled manner. In this paper we 
will consider only two-dimensional motion (p z = z = Z\ = z 2 = 0, say) with the further 
specializations X\ = 0, x 2 = (x p ,0), and Gm\ = 1 ^> Gm 2 , so that the motion is nearly- 
Keplerian except very near m 2 . There are three classes of bound motion (unbound orbits 
will not be considered). First, the particle can be tightly bound to either m\ or m 2 , never 
closely approaching the other mass; this is the limit in which the motion is (quantitatively) 
similar to that in the Stark problem. In the second type of motion, the test body is confined 
to the annulus between two confocal ellipses (with foci at the positions of mi and m 2 ) 
and hence maintains a finite distance from both masses at all times. In the third type of 
motion, the inner bounding ellipse disappears and the particle eventually approaches each 
mass arbitrarily closely and with encounter eccentricities arbitrarily close to unity. Examples 
of the latter two motions are shown in Figure || 



We shall focus attention on the integration of initial conditions lying near the critical 
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Fig. 8. — Two examples of bounded motion in the two fixed point problem (§ |5.1| ). In 
the top panel, the motion is confined to lie between two ellipses (with foci at the indicated 
positions of the fixed masses) and hence maintains a finite distance from both masses. In the 
bottom panel, the inner ellipse no longer exists and arbitrarily close (and radial) encounters 
with both masses are possible. 
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line between motion of the second and third types; on the critical curve, the inner ellipse 
degenerates to the line connecting mi and m 2 and the limiting motion consists of an ever- 
tightening spiral which converges on this line. Convenient, precisely critical initial conditions 
are Xq = (1, 0), p = (0, 1), and x p = —1 (m2 remains arbitrary); slight variations of x p from 
this critical value then allow exquisite control over the severity of close encounters. 



5.1.2. Numerical Integrator Performance 

Numerical results are shown in Figures P[|TT|; the relative execution times for each simu- 
lation are given in Table [TJ Each figure contains six panels, one for each integration method 
tested: the original Wisdom-Holman method (labeled 'WH'), the regularized WH mapping 
('RWH', § j37T|) , the original potential-splitting method ('PS'), our modified PS algorithm 



MPS', § |3.2|) , the time-reversible Stark scheme ('TRS', § |4.1|) , and the regularized, sym- 



plectic Stark method ('RSS', § |4.2|) . The Stark stepper for the TRS and RSS integrators 
implemented the analytic solution in the case of bound motion; for unbound motion, a highly 
accurate (machine level truncation error) Bulirsch-Stoer routine was used to numerically in- 
tegrate the K-S regularized equations of motion. Each panel plots the relative errors in 
the energy E (solid curve) and separation constant a (dotted line) over the course of an 
integration lasting 10 3 or 10 4 orbital periods, where t or b is the period of the unperturbed 
Kepler orbit. (Since in most cases the errors in E and a are very similar, the two curves 
are often difficult to distinguish.) All integrations used At = i0 -3 t O rb, Gm 2 = 0.01, and 
the aforementioned near-critical initial conditions; only the value of x p was changed between 
figures. 

Figure |9] plots the results for x p = — 1.5. In this case there are never close encounters 
with m2 and the motion is similar to that in the Stark problem; in particular, the orbit 
periodically becomes highly eccentric. The qualitative results in this case are quite simple — 
all regularized schemes (RWH, MPS, RSS) behaved well and showed no signs of instability, 
whereas all unregularized ones proved to be unstable at high eccentricities. In addition, 
although the TRS scheme performs well initially (which is not surprising, since the motion 
is Stark-like) its errors grow secularly and by the end of the simulation are of order unity; 
this is disappointing considering that the method is explicitly time-symmetric, and further 
supports the 'non-symplectic' label previously placed on it (see § PO). 



Figure [10| shows the results for x p = —1.02. In this case there are frequent close 
encounters with both masses and there is no limit to how close and radial they can be (cf. 
Figure||, lower panel). As expected, all the single-timestep integrators (WH, RWH, TRS, 
RSS) quickly falter, unable to cope with the close encounters. Although initially stable, the 
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Fig. 9. — Comparative integrator performance for the two fixed point problem with x p = 
—1.5 and Gm<i = 0.01 (see § 5.1.2]) . In this case there are no close encounters with the mass 
m 2 and the orbital motion is similar to that in the Stark problem (cf. Figure [I]). Each panel 
displays the relative errors in energy E (solid curve) and separation constant a (dotted curve) 
for a specific integration scheme. All regularized methods (RWH, MPS, RSS) are stable and 
well-behaved; all others go unstable when the orbit becomes nearly radial. The TRS scheme 
also exhibits linear error growth. 
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Fig. 10. — Similar to Figure p] but for x F 



-1.02. Here the orbital motion is like that 



in the lower panel of Figure |^ and the test particle undergoes arbitrarily close (and radial) 
encounters with both masses. As expected, all the single-timestep integrators (WH, RWH, 
TRS, RSS) quickly falter since they cannot handle close encounters. Although initially 
stable, the PS integrator (here including regularization) is occasionally overwhelmed by the 
extremely close encounters present in the problem. Only the MPS algorithm, incorporating 
both regularization and force-center switching, performs satisfactorily. 
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Fig. 11. — Similar to Figure |9] but for x p = —0.95. The orbital motion now corresponds 
to the top panel of Figure ||: both masses are approached, but encounter distance and 
eccentricity both have finite bounds (b ^ 0.05 and e <; 0.95, respectively). The chosen 
timestep, At = 10" 3 t or b, was sufficient to resolve all close encounters and hence all integrators 
(except the TRS method) perform well. As in Figure |S], the TRS scheme displays linear error 
growth. 
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PS integrator is occasionally overwhelmed by the extremely close encounters present in the 
problem; further investigation revealed that the discrete jumps occurred during encounters 
with impact parameters b <^ 10~ 3 . (In this instance the PS integrator was regularized 
since otherwise it would have been unstable even in the absence of close encounters.) Only 
the MPS algorithm, which also included force-center switching, performs satisfactorily. As 
evidenced by the formation of a secular trend late in the test, however, even this algorithm 
has its limits. In particular we have not found a practical way to regularize about the 
perturber while the force-center switch is active, and hence the scheme is unstable whenever 
the encounter eccentricity exceeds some limit. The growing error may also be partly due to 
non-symplectic behavior (or 'ringing') introduced during the switching process, but as both 
effects occur only during encounters it is difficult to distinguish between the two (at least in 
these simulations). 



The results for x p = —0.95 are displayed in Figure [IT]. Here the qualitative orbital 
motion corresponds to the top panel of Figure § although both masses are approached, all 
encounters have impact parameters b ^ 0.05 and local eccentricities e <^ 0.95. This implies 
that the chosen timestep (At = 10~ 3 t or b) is just adequate to clearly resolve the encounters 
with both mi and m 2 . In this case, therefore, nearly every integrator performs quite well. 
The exception is the TRS method, which is again plagued by a strong linear error growth. 
The PS and MPS routines do particularly well because of the effectively smaller timestep 
used near m 2 (the approaches to m 2 are close enough that the algorithms subdivide the 
timestep several times each encounter), but the relative cost in execution time (see Table [[]) 
is commensurate — the WH and RWH schemes would have produced comparable or superior 
results given the same amount of CPU time. 

The timing results listed in Table |l]are straightforward to interpret. The RWH method 
is the most efficient due to the speed of the regularized Kepler stepper, the force calculation 
being essentially trivial here. The PS and MPS methods show little overhead cost in the 
absence of close encounters, but slow down substantially when close approaches occur fre- 
quently. The rather obvious conclusion is that multiple-timestep integrators (and other 'en- 
counter codes') should only be used when the detailed encounter dynamics are important — if 
a moderately softened perturbing potential is acceptable, for example, integration using a 
constant timestep scheme would be both stable and more efficient. The Stark-based schemes 
are at least an order of magnitude slower than the corresponding Kepler mappings because of 
the high cost of taking Stark steps instead of Kepler steps, yet in this case were no more accu- 
rate. We conclude that such methods are uncompetitive unless the perturbation potential is 
extremely Stark-like; in particular, for point mass perturbers — which only appear Stark-like 
at large distances and to lowest order — we consider them to be of marginal interest. 
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5.2. Galactic Dynamics Test Problem 

The results of a more generic test problem are shown in Figure [12] (relative timings can 
again be found in Table [I]). The problem consisted of integrating test particle motion in 
the perturbing field of 100 fixed points of mass 10 _3 M (where M is the mass of the central 
object, held fixed at the origin). The positions of the masses were drawn randomly from a 
spherically symmetric distribution with a radial density profile oc r™ 2 ; this is similar to the 
mass distribution seen in several galactic nuclei believed to contain massive black holes, and 
closely resembles the configurations used by Rauch & Ingalls (1997). (We do not allow the 
perturbers to orbit M because energy would not be conserved in this case.) In addition, 
for the WH, RWH, TRS, and RSS integrators (i.e., the constant timestep methods) the 
perturbing potentials were slightly softened to limit the severity of close encounters and 
allow a more realistic comparison with the multiple-timestep (PS and MPS) routines to be 
made. The step size in all cases was At = 10 -3 t or t,. Finally, we note that the MPS integrator 
used in these simulations did not include force-center switching, which would have required 
substantial (though conceptually straightforward) modifications to the original code due to 
the multitude of perturbers involved. 

As before each panel in the figure plots the relative energy error of a particular integrator. 
The solid line represents the results for an orbit with low initial eccentricity (eo = 0.5), 
and the dotted line for one with high eccentricity (eo = 0.99); the initial conditions were 
otherwise identical and were the same for every integrator. The results are completely in line 
with previous ones and show nothing unexpected. Only the regularized methods are stable 
at high eccentricities, the PS and MPS algorithms appear highly successful at resolving the 
many (unsoftened) close encounters that occur, and once more the error in the TRS scheme 
is dominated by a secular drift. 

6. Discussion 

We have shown that the WH mapping is generically unstable when applied to eccentric, 
nearly-Keplerian orbits whenever the step size is not small enough to resolve periapse. This 
'radial orbit instability' is fully explainable in terms of the overlap of step size resonances 
and has a simple geometric manifestation in the case of the Stark problem. Our investigation 
indicates that the islands of stability found in the latter problem do not exist in the more 
general cases we have examined; the instability therefore appears to be unavoidable in typical 
situations, unless one employs the brute-force approach of decreasing the timestep by the 
requisite amount. However, besides the fact that this is an extremely inefficient solution — it 
reduces the mapping to a very costly direct integration scheme — we have shown that an 
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Fig. 12. — Results for the galactic dynamics test problem (§ |5T2| ; cf. Figures |9"14TTD. The 
solid curve shows the energy error for an orbit with low initial eccentricity (eo = 0.5); the 
dotted line corresponds to eo = 0.99 (the rest of the configuration remaining unchanged). The 
performance of each integrator is in line with expectations; in particular, only the regularized 
methods remain stable when e ~ 1, and the TRS method continues to be dominated by a 
secular trend. 
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elegant solution to the problem is already available: the regularization approach of Mikkola 
(1997). In every case examined, not only was the regularized WH mapping immune to 
the radial orbit instability, in many cases it was also more efficient. We enthusiastically 
recommend its use whenever close encounters with perturbers are not of concern. 

We remind the reader that our investigation has not cast into doubt all previous studies 
that have used the WH integrator and its variants. In nearly every case care has been taken 
to use a small enough step size such that perihelion passage would be adequately sampled. 
We note only one area where particular caution should be exercised. One of the features of 
the long-term dynamics in mean motion resonances and secular resonance is that very high 
eccentricities can be developed. These eccentricities are often large enough that physical 
collision with the sun is a common outcome in studies of meteorite delivery from the main 
asteroid belt and the long-term dynamics of ecliptic comets ( |Gladman et al. 1997| ; [MorbideTT 



Moons 1995; Lcvison fc Duncan 1997). In those cases it is unlikely that the step size used 



was small enough to resolve the perihelion passage. Although these researchers checked their 
results for step size dependence and reported no numerical artifacts, we suggest that further 
examination of those cases would be prudent. 

We have demonstrated that the potential-splitting method of Lee et al. (1997) can 
be regularized to produce an algorithm that is robust in the face of both close encounters 
and highly eccentric orbits. We have also shown that force-center switching during excep- 
tionally close encounters can be cleanly incorporated into the method and can substantially 
enhance the stability of the algorithm without noticeably affecting its desirable symplectic 
qualities. We have not, however, found a practical way to regularize around the perturber 
while the switch is in effect; the stability of this approach during highly eccentric encounters 
is correspondingly questionable. 

Our examination of Stark-based integrators indicates that they, too, are subject to 
the radial orbit instability unless regularized, although it tends to be less severe since the 
Stark approximation becomes systematically better near the origin. Unless the perturbing 
potential is very well represented by a Stark potential, they also appear uncompetitive in 
terms of efficiency — by over an order of magnitude — due to the cost of Stark steps relative to 
that of Kepler steps. Among the Stark-based methods, the regularized, symplectic method 
(§ f4.2| ) consistently outperformed the time- reversible method (§ [4.1| ), in part because of the 
linearly growing energy error exhibited by the latter. Our conclusion is that Stark-based 
schemes are of marginal utility in the integration of N-body systems. 

It is clear that integrators based on a two fixed point (TFP) splitting instead of a Stark or 
Kepler splitting are also possible; they can be constructed in the same manner as the Stark- 
based methods were. Such methods could be useful whenever two bodies strongly dominate 
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the mass in the system (e.g., asteroid motion in the Sun- Jupiter system). As for Stark 
motion, however, the relative expense of advancing the TFP Hamiltonian is a significant 
handicap, and the circumstances in which its use is justified remain unclear. On the other 
hand, since Stark motion is a subset of TFP motion it is not unlikely that methods based on 
the latter splitting will generically outperform those of the former type, since their analytic 
solutions are of similar complexity. It would be interesting to investigate this possibility in 
greater detail. 

Although we have confined attention to the perturbed two-body problem, the techniques 
employed in this paper are also applicable to general hierarchical iV-body systems. In partic- 
ular, we believe that regularization of the iV-body version of the potential-splitting method 
( puncan, Levison, fc Lee 1997| ) is likely to cure the instability at high eccentricities noted 
by the authors. In principle force-center switching of the kind described in § FO can also 



be done, but we have not studied this possibility in detail. In any event, we have found the 
combination of regularization and potential-splitting to be a powerful one, and to produce 
a remarkably versatile symplectic method for the integration of nearly-Keplerian systems. 



We thank Doug Hamilton, Norm Murray, Scott Tremaine, and Man Hoi Lee for illumi- 
nating discussions. 
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Table 1. Relative Efficiency of Selected Integration Simulations 



Integrator 



Fig. | 



Fig. [TO 



Fig. P 



Fig. IT2 



WH 
RWH 
PS 
MPS 
TRS 
RSS 



1.00 
0.75 
1.20 
1.25 
12.8 
13.7 



I. 00 
0.73 
14.5 
20.3 
15.0 

II. 9 



1.00 
0.74 
6.40 
8.50 
20.0 
13.8 



1.00 

I. 40 
17.5 
25.5 
12.8 

II. 1 



